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Abstract. For some physics students, the concept of a particle travelling faster 
than the speed of light holds endless fascination, and Cerenkov radiation is a visible 
consequence of a charged particle travelling through a medium at locally superluminal 
velocities. The Heaviside-Feynman equations for calculating the magnetic and 
electric fields of a moving charge have been known for many decades, but it is 
only recently that the computing power to plot the fields of such a particle has 
become readily available for student use. This article investigates and illustrates the 
calculation of Maxwell's D field in homogeneous isotropic media for arbitrary, including 
superluminal, constant velocity, and uses the results as a basis for discussing energy 
transfer in the electromagnetic field. 
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1. Introduction 

For any medium, the speed of light is also the speed of propagation of electromagnetic 
waves within that medium — after all, light is simply an electromagnetic wave. But 
what is the relationship between an electromagnetic wave and the electromagnetic field 
of a particle? 

In answer to this, we consider the field of a stationary charge in free space. The 
field of the charge q extends to infinity and at distance r is simply given by 



If, however, we consider the field to be made up of an imaginary flux originating at 
the particle itself and travelling radially outwards with velocity c, we see that the flux 
present at distance r at time t will have originated at the particle at a time 



When the particle is stationary, it is unnecessary to make this distinction as the point 
from which the field was emitted is the same as the point currently occupied by the 
particle. When the particle is in motion, however, it becomes necessary to identify both 
when and where the particle may have been located so as to give rise to the field at a 
particular point. 

Let the particle follow an arbitrary path x(t). We may write the four- dimensional 
co-ordinates of any point on this path in the form (x(t),t). We now wish to determine 
the field at point (x 1 , t'), which may or may not lie upon the path. In order for the field 
to have propagated here at speed c, it must have originated at a co-ordinate (x,t) such 
that 



A value of (x(t),t) satisfying (3) is referred to as the retarded position of the particle 
generating the field at (x',t'). For certain trajectories, such as superluminal motion, 
there may be multiple retarded positions contributing to the field at a given location. 
These are combined using the principle of linear superposition. 

So how do we determine the fields arising from the particle at the retarded position? 
The particle is in motion, and may even be accelerating. Obviously the Coulomb field is 
not appropriate; the moving particle constitutes a current as well as an electric charge. 
Writing R = x' — x and R = \R\, the equation for the relevant magnetic field 




(2) 



x — x \ — 



c(t'-t). 



(3) 





x R 



Visualisation of Cerenkov Radiation and the Fields of a Moving Charge 
Robert N. C. Pfeifer and Timo A. Nieminen 
European Journal of Physics 27 (2006) 521-529 



3 



is attributed to Oliver Heaviside [1, p. 436], and the explicit equation for the electric 
field 
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although also originally developed by Heaviside [1, p. 437], is usually attributed to 
Richard Feynman [2, p. 11-21- 1]. They are related by 

B = R x E. (6) 

The corresponding scalar and vector potentials are known as the Lienard-Wiechert 
potentials 

' 1 q 



<j>(x',t) - 
A(x',t) 
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(7) 



(8) 



Derivations of the above expressions may be found in advanced electromagnetics 
textbooks such as Jackson [3]. 

Now let us extend our considerations to media other than vacuum. In material 
media, the values of e, fi and refractive index n will generally differ from those in 
vacuum, and the medium may also be anisotropic. Furthermore, when a charge is in 
relativistic motion the medium will appear to contract along the direction of travel, with 
consequences for e and fi. Thus even a medium which is isotropic at rest will appear 
anisotropic to a moving charge. 

Our calculations shall be performed in the rest frame of an isotropic medium. 
However, by working with Maxwell's D and H fields we take an approach more readily 
extensible to the cases of moving and anisotropic media. 

Because 



and 



D = sE 



B = fiH, 



(9) 



(10) 



e and fi are eliminated from our equations and only the value of the refractive index 
affects our results. This alters the speed of propagation of field modulations within the 
medium, 

c 



c = 



n 



(11) 
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and thus affects the solutions for the retarded positions. The equations for the D and 
H fields then become 

d 








R 2 ) 






/ ret 




xx R' 
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^ ret 




ret 
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(12) 



(13) 



Note that the c associated with each is unchanged as it arises not from the field 
propagation time, but as a constant relating our units of measurement in temporal and 
spatial dimensions. This is unaffected by a change of material medium. 



2. Calculation of Retarded co-ordinates 



Consider a particle travelling in a straight line parallel to the y-axis at constant velocity 
v. We choose our co-ordinate system such that the direction of travel is along the x-axis 
and wish to calculate the D field magnitudes in the xy plane. If the particle is at y 
co-ordinate y at time t then its equation of motion is given by 

y = y + v(t-t ). (14) 

Suppose we wish to calculate the D field at spatial co-ordinate {x' , y', 0) at time t! . 
We denote particle co-ordinates with unprimed characters and field co-ordinates with 
primed characters. As can be seen from figure 1, the following relationships hold: 

R _ Rn 

d c 



(f - 1) 



(15) 



R 2 



x' 2 + (y'-y) 2 (16) 

Combining (14), (15) and (16), and choosing c = 1 and t = for clarity, we obtain a 
quadratic in y: 

y 2 (n 2 v 2 - 1) + y(2y + 2t'v - 2y'n 2 v 2 ) 

+ ( x ' 2 n 2 v 2 + y' 2 n 2 v 2 -t' 2 v 2 -y 2 -2t'y v)=0. { ' 

In accordance with convention we discard solutions corresponding to advance potentials 
(t > t'), using (14) to identify the time t corresponding to a given source solution y. For 
subluminal particles in a single medium, one of the solutions will be advanced and the 
other retarded. For superluminal particles, solutions will either be both advanced, both 
retarded, or both imaginary, indicating that the field has not yet reached this region. 

Having identified the space-time co-ordinates of our retarded source(s), we can 
now calculate the field at (x',y',0,t') using (12). When this process is repeated for 
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multiple sets of co-ordinates (x', y', 0, t'), a plot of the field may be built up. Figures 2- 
4 demonstrate the results which may be obtained. 

3. Results 

To represent the field of the particle, we have plotted contours of equal D field 
magnitude. This should not be confused with the plotting of electric or magnetic field 
lines. The reader may be familiar with the representation of a Coulomb field undergoing 
a Lorentz boost shown in figure 5. Why, then, is the field plotted in figures 2 and 3 not 
similarly symmetrical? 

When the charge is placed in motion, it constitutes a finite current element and will 
therefore generate a magnetic field. The charge is moving, and hence distances to this 
current element will vary with time. As they do, the magnetic field will also vary with 
time, inducing a further electric field. This field is opposed to the existing compressed 
Coulomb field in the regions marked 'a' in figure 5, and complements it in the regions 
marked 'b', giving rise to the observed distortion in the field contours. The magnitude of 
this distortion is dependent on electromagnetic induction in accordance with Maxwell's 



5 



Visualisation of Cerenkov Radiation and the Fields of a Moving Charge 
Robert N. C. Pfeifer and Timo A. Nieminen 
European Journal of Physics 27 (2006) 521-529 



6 




X 



Figure 2. Field of a charged particle in a homogeneous medium; 
n = 2, v = 0.45 c 

equations, and hence depends on the speed of the particle relative to the speed of light 
in vacuo. Hence less distortion is noted in figure 2 than in figure 3, despite the particle's 
speed being 0.9 times the speed of light in the local medium in each case. 

In figure 4, the particle is now travelling superluminally and as a result continually 
overtakes the leading edge of its propagating field. Because of this, a shock front is built 
up. It is this which is perceived as the Cerenkov radiation (see section 4). An analogy 
which is often employed to illustrate this phenomenon is to imagine that as it travels, 
the particle emits repeated pulses of electromagnetic radiation, expanding in shells as 
shown in figure 6. While a useful aid to visualisation, this analogy breaks down when 
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Figure 3. Field of a charged particle in a homogeneous medium; 
n = 1, v = 0.9c 

applied to subluminal particles and regions lying within the Cerenkov cone. A more 
complete explanation could be developed in which the field constitutes the emission of 
virtual photons in analogy with quantum field theory, but lies outside the scope of this 
article. 

It is also interesting to see how these results relate to the time-reversal symmetry 
of the Maxwell equations. Essentially, time reversal interchanges the advanced and 
retarded potentials, and hence appropriately reverses the direction of the Cerenkov cone, 
as would be expected if the direction of motion of the charge were reversed. This can 
be contrasted with the Lorentz-contracted field lines of figure 5, which are intrinsically 
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Figure 4. Field of a charged particle in a homogeneous medium; 
n = 2, v = 0.9 c 



symmetric, and for which time reversal simply reverses the direction of travel of the 
particle leaving the field lines unchanged. 

4. Discussion 

4-1. Electromagnetic Radiation 

As is clearly shown in figure 7, the shock cone constitutes a narrow region of 
comparatively powerful electric and magnetic fields, giving rise to local concentrations 
in field energy. These regions of concentrated field propagate outward with time. Of 
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Direction of motion 




Figure 5. Field line representation of a coulomb field undergoing a Lorentz boost. 
The particle is travelling left-to-right so regions marked 'a' lie ahead of the moving 
particle with respect to the direction of motion, and those marked 'b' lie behind it. 



Figure 6. Formation of a shock cone from repeated emission of spherical wavefronts 

course, a freely propagating energy- carrying wave in the electric and magnetic fields is 
what we know as electromagnetic radiation, in this case visible light. 

In allowing the particle to continue to pursue a constant velocity trajectory in our 
calculations, we have neglected the effects of this radiative energy loss. 

4-2. Freedom to Propagate 

What constitutes a freely propagating electromagnetic wave? In the above section we 
identified the field surge of the Cerenkov cone with electromagnetic radiation. But the 
method used to generate these images makes no allowance for free packets of fields 
propagating through space — the fields involved all originate directly from the moving 
charge. Is it therefore appropriate to think of this wave as free in the same sense as we 
think of photons being free? 




9 



Visualisation of Cerenkov Radiation and the Fields of a Moving Charge 
Robert N. C. Pfeifer and Timo A. Nieminen 
European Journal of Physics 27 (2006) 521-529 



10 




Figure 7. '.D' field of a superluminal particle (see also figure 4) 

The answer is yes. Although in particulate models photons are considered 
as independent entities, and likewise in classical electromagnetics we often consider 
sourceless, freely propagating plane waves, in practice there exists a charge at the end of 
every photon or electromagnetic wave, of whose retarded fields it is in fact an extension. 
This charge may be accelerating, jumping between atomic orbitals, or travelling through 
an optically dense medium as seen here. Our 'free wave' is indeed free, in that its nature 
and behaviour are unaffected by any subsequent actions of the originating charge: That 
modulation in the local electromagnetic field will continue to propagate out indefinitely 
at the local speed of light, even if the originating particle is subsequently somehow 
destroyed. 

4-3. Acoustic Shockwaves 

As we have seen in section 3, Cerenkov radiation arises due to the formation of a 
shockwave in a particle's electric and magnetic fields. The analogous problem in 
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Figure 8. Supersonic gas flow over a stationary cone at Mach 4, visualised using the 
background oriented schlieren technique. 

acoustics is the formation of Shockwaves due to supersonic gas flows over material bodies, 
which continues to be of great importance in aeronautical engineering. In 1886 the first 
photographs of the bow shock of a supersonic projectile were created by the collaboration 
of Mach, Salcher and Riegler [4], utilising a technique derived by Toepler in 1864, known 
as the schlieren method. Related techniques continue to be in use to this day (figure 8). 

Once again, the shock cone may be considered to be built up by superposition of 
consecutive spherical wavefronts emitted by the source as it travels (figure 6). In this 
case, the wavefronts are pressure waves within the surrounding medium. Diagrams 
showing the construction of the shock cone in this manner were first published by 
Christian Doppler in 1842 [5], though it is after Mach that the cone is usually named, 
in recognition of his later experimental work. 

Introduction of factors such as viscosity and turbulent flow lead to additional 
behaviours not discussed in this article, and cause the study of fluid dynamics and 
acoustic Shockwaves to be a complex and fascinating field. 



11 



Visualisation of Cerenkov Radiation and the Fields of a Moving Charge 12 
Robert N. C. Pfeifer and Timo A. Nieminen 
European Journal of Physics 27 (2006) 521-529 

5. Supplementary Material 

The electronic version of this article is accompanied by the MATLAB program 
Cherenkov . m, which was used to generate the plots accompanying this article. Usage 
instructions may be viewed by typing 'help Cherenkov .m'. 

6. Conclusion 

This paper has aimed to illustrate the behaviours of the fields of a moving charge in 
an optically dense medium, including Cerenkov radiation. These behaviours are readily 
simulated on a modern desktop computer, and demonstrate how the radiation of the 
Cerenkov cone arises naturally from the fields of a superluminal charge. The concept of 
a 'free' photon is discussed, in relation to its origin in the retarded field of an electric 
charge. It is explained how the existence of the photon arises as a result of the motion 
of the charge, but that the subsequent behaviours of the charge and the photon, or field 
wave packet, are independent. Finally, an analogy is drawn between the formation of 
the Cerenkov cone and the formation of the bow shock of a supersonic projectile, a topic 
of vital importance and ongoing research in aeronautical engineering. 

It is hoped that this exhibition of the interesting phenomenon of Cerenkov radiation 
may stimulate the student to further self-guided learning, whether by developing upon 
the theme of this article (for example, by simulating the fields of a particle pursuing an 
arbitrary path, or adjacent to a medium of differing refractive index) or by investigating 
other specific radiative phenomena. For example, the Heaviside-Feynman equations 
may be used to calculate the fields of an accelerating charge as found within a 
radiating dipole antenna, or a synchrotron. An advanced student pursuing an interest 
in astrophysics may wish to investigate the radiative consequences of a straight particle 
path in curved space-time. There are many more interesting possibilities to explore. 
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Software 

Cherenkov.m 

function Cherenkov (ymax , yf in , tf in , v , ref index , t , plotopt ions) 
7o Cherenkov (ymax , yf in , tf in , v , ref index , t , plot opt ions ) 

7, Calculates Maxwell's D field for a charged particle travelling through a 
7 homogenous medium. 

7« ymax sets the size of the viewing area. 

°/« Particle travels along line y=ymax/2, reaching yfin at time tfin. 
7« Velocity v is a decimal multiplier of c. 
7, Plot is generated for the field at time t. 

7. 

'/, plotoptions: 
7. 1: 3D 

7o 2: coloured surface 
7. 4: contour (slow) 

7o Options may be combined using addition. 
7, Default value: 7 (all). 

7. 

7o Examples : 

7. Try Cherenkov (150, 120, 100, 0.8, 2, 100, 3) 
7. and Cherenkov (150, 120, 100, 0.8,1, 100) . 

if (nargin==6) 

plotoptions=7; 

end 

if (nargin<6) 

'Insufficient arguments supplied. Please type help CherenkovP for more information.' 
return 

end 

warning off MATLAB : divideByZero 

particlepos=ymax/2 ; °/„ path goes along this x line 
advret=+l; 7. +1 or -1 for advance or retarded potentials 
yzero=yf in-v*tf in; 7 particle position at time t=0 

7o Need to tabulate unit R vector, R magnitude and lookback time 
7o for all 4-co-ords wrt particle track. 

R=zeros(ymax,ymax,3,3) ; 7o x,y,t, values 

Rvec=zeros(ymax,ymax,3,2,2) ; 7. x & y unit vectors - 2 sets of solutions. 
7o First three indices indicate location at which solutions apply. Fourth 
7» index indicates x and y components of solution. Fifth index indicates 
7o first and second valid solutions. Y solution space is also used as 
7, temporary workspace . 
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7 Fill R(,l) with x positions, R(,2) with y positions, R(,3) with z 

°/« positions, R(,4) with t positions; i.e. make R the position vector array. 

R(l :ymax, 1 , 1 , 1)=1 : ymax; 

R(l,l:ymax,l,2)=(l:ymax) ' ; 

for q=2:ymax 

R(: ,q,l,l)=R(: ,1,1,1); 

R(q,:,l,2)=R(l,:,l,2); 

end 

R(:,:,2,1:2)=R(:,:,1,1:2); 
R(:,:,3,1:2)=R(:,:,1,1:2); 
R(: , : ,l,3)=t-l; 
R(: , : ,2,3)=t; 
R(: , : ,3,3)=t+l; 

7« Put working value 1 into Rvec(,l,l) 

Rvec(: , : , : ,1,1)= (yzero + R(:,:,:,3)*v - v . ~2 . *R( : , : , : ,2) . *ref index. "2) ; 
Rvec(: , : , : ,1,1)= Rvec(: , : , : ,1,1) ./ (v . ~2 . *ref index . ~2 - 1) ; 

7 Put working value 2 into Rvec(,l,2) 

Rvec(: , : , : ,1,2)= ref index . ~2 .* (R( :,:,:, 2) . ~2+(particlepos-R( :,:,:, 1) ). ~2) - R( : , : , : ,3) . ~2; 
Rvec( : , : , : , 1,2)= v. ~2.*Rvec( : , : , : , 1 ,2) - yzero. ~2 - 2*yzero*v*R( : , : , : ,3) ; 
Rvec(: , : , : ,1,2)= Rvec(: , : , : ,1,2) ./ (v . ~2 . *ref index . ~2 - 1) ; 

7o Put source Y solution 1 into Rvec(,2,l) 

Rvec( : , : , : ,2, 1)= -Rvec( :,:,:, 1 , 1) + sqrt( Rvec( :,:,:, 1 , 1) . "2 - Rvec (:,:,:, 1 , 2) ); 
7o Put source Y solution 2 into Rvec (,2,2) 

Rvec( : , : , : ,2,2)= -Rvec( :,:,:, 1 , 1) - sqrt( Rvec( :,:,:, 1 , 1) . "2 - Rvec (:,:,:, 1 , 2) ); 

7 Put retarded source time into sourcet 
sourcet(: , : , : ,1, :)=(Rvec(: , : , : ,2, :) -yzero) /v; 

7o Change from source absolute Y to relative Y 
Rvec(: , : , : ,2,1) = R(:,:,:,2) - Rvec( : , : , : ,2, 1) ; 
Rvec(: , : , : ,2,2) = R(:,:,:,2) - Rvec( : , : , : ,2,2) ; 
7« e.g. +20 means you are 20 ahead of your source point 

7. Store retarded source X displacement in Rvec(,l,:) 
Rvec(: , : , : ,1,1) = R( :,:,:, 1) -particlepos ; 
Rvec(: , : , : ,1,2) = R( :,:,:, 1) -particlepos ; 

7.sourcey=Rvec (:,:,:, 1,:); 

7o R is now going to be displacement magnitude to retarded source: 
clear R 

R=zeros(ymax,ymax,3,2,2) ; 
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°/« Put distance from retarded source into R(,l) and R(,2): 
R(: , : , : ,1, : )=sqrt (sum(Rvec . ~2 ,4) ) ; 
R (:,:,:, 2 ,:) =R (:,:,:, 1 ,:) ; 

°/« Divide Rvec by R to make it a unit vector 
Rvec=Rvec . /R; 

7, Create array for Rvec on R~2 term 
RonRsq=Rvec . /R. /R; 

'R matrices done' 

Df ield=RonRsq+ddt (RonRsq) . *R+ddt (ddt (Rvec) ) ; 

°/ Eliminate source distances of zero (which create NaN entries which are 
°/ most troublesome) 

for xx=l:ymax 

for yy=l:ymax 
for sol=l:2 

for coord=l:2 

if isnan(Df ield(xx,yy,2,sol,coord)) 
Df ield(xx,yy,2,sol,coord)=0; 

end 

end 

end 

end 

end 

'D field done' 

7, Reduce array of vector D field to array of scalar magnitude of D field 
Dfield=sqrt(sum(Dfield.~2,4)) ; 

7« Eliminate invalid sources: 

if (advret==l) 

°/ Eliminate sources in past : 
Df ield=Df ield. *(sourcet<t) ; 

else 

°/, Eliminate sources in future: 
Df ield=Df ield.*(sourcet>t) ; 

end 

7« Eliminate contributions due to imaginary particle times/positions (! ) 
Df ield=Df ield. *(imag(sourcet)==0) ; 
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'Unwanted sources eliminated' 

7 Combine valid sources by linear superposition 
Df ield=sum(Df ield,5) ; 

7. Orient graph correctly for plotting 
DfieldC: , : ,2, : )=Df ield( : , : ,2, :) ' ; 

if (bitand(plotoptions ,4) ) °L Contour plot 
'Creating contour plot. . . ' 
figure (3) 

v= [0:0. 0001: 0.01] ; 

contour (real (DfieldC : , : ,2 , : ) ) , v) ; 

colormapC [0 0]); 

xlabel('x') 

ylabel('y') 

end 

if (bitandCplotoptions , 1) ) 7, 3D plot 
'Creating surface plot. . . ' 
figure (1) 

surf c (real (DfieldC : , : ,2, :))) ; 

caxis( [0 0.02] ) ; 

zlimCCO 0.2]); 

xlabelC'x') 

ylabelC'y') 

zlabelC'D') 

end 

if CbitandCplotoptions ,2) ) % Flat colour plot 
'Creating colour plot. . . ' 
figure C2) 

pcolorCrealCDfieldC: , : ,2, :))) ; 

shading flat 

caxisC [0 0.02] ) ; 

xlabelC'x') 

ylabelC'y') 

end 



ddt.m 

function B=ddtCA) 

7, Takes derivative wrt time of field 

[x,y,t,s,d]=sizeCA) ; 
BiggerA=zeros Cx , y , t+2 , s , d) ; 
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B=zeros(x,y,t,s,d) ; 

BiggerA(l:x,l:y,2:t+l, : , : )=A(1 :x, 1 :y, 1 :t, : , :) ; 

"/.Linearly extrapolate beyond boundaries of original A to first order 

BiggerA(l:x,l:y,l, : , : )=A(1 :x, 1 :y , 1 , : , : ) . *2-A(l :x, 1 :y ,2, : , :) ; 
BiggerA(l:x,l:y,t+2, : , : )=A(1 :x, 1 :y,t, : , : ) . *2-A(l :x, 1 :y ,t-l , : , :) ; 

"/.Calculate ddt's 

B(l:x,l:y,l:t, : , :) = ( BiggerA(l :x, 1 :y,3:t+2, : , : ) -BiggerA(l : x, 1 :y , 1 : t , : , :) ) ./2; 
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